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The relation between the Lattice Boltzmann Method, which has re- 
O . cently become popular, and the Kinetic Schemes, which are routinely used 

in Computational Fluid Dynamics, is explored. A new discrete velocity 
model for the numerical solution of Navier-Stokes equations for incom- 
pressible fluid flow is presented by combining both the approaches. The 
new scheme can be interpreted as a pseudo-compressibility method and, 
for a particular choice of parameters, this interpretation carries over to the 
Lattice Boltzmann Method. 



1 Introduction 

In the last decade, Lattice Boltzmann Method has emerged as a potential alterna- 
tive to other Computational Fluid Dynamics techniques in simulating fluid flows 
numerically. The Lattice Boltzmann Method (LBM) was first introduced by Mc- 
Namara and Zanetti |lj to overcome the drawbacks of the Lattice Gas Cellular 



f 



Automata (LGCA), which resulted from attempts to obtain macroscopic fluid 
flow simulations from the simplest possible microscopic description using discrete 
velocity models. See the references []2], |3[] and Q for reviews of the LGCA, refer- 
ences and for reviews of the Lattice Boltzmann method and |7J for a review 
on discrete velocity models. 

Some authors noted the closeness of the Lattice Boltzmann Method to the 
Kinetic Schemes (see [BJ, M), which have been routinely used in CFD simulations 
(see reference |L0| for a review of Kinetic Schemes). Both methods use the Boltz- 
mann equation of Kinetic Theory as the starting point, but are aimed at solving 
the macroscopic equations of fluid flow. This approach exploits the fact that the 
Boltzmann evolution is essentially equivalent to Euler or Navier Stokes evolution 
if the state is in or close to local thermodynamic equilibrium. While most of 
the Kinetic Schemes were developed for the solution of compressible equations, 
the Lattice Boltzmann Method operates in the incompressible limit. However, 
as we will show in Section |2|, the two methods even coincide for a particular 
parameter constellation. This implies that the observations which are valid for 
Kinetic Schemes also have a direct consequence on LBM. In view of these re- 
marks it is surprising that the close relation between the two methods is not fully 
appreciated. Our intention is to stress the remarkable coincidence. 

The first Kinetic Scheme was introduced more than two decades back, by 
Sanders and Prendergast [11]. It is popularly known as the Beam scheme, and, 



incidentally, is also a discrete velocity model for simulating Euler equations. A 
few years later, an approach to construct Kinetic Schemes for general hyperbolic 
systems of conservation laws was described by Harten, Lax and van Leer [|T^] . 
Many Kinetic Schemes for the compressible Euler system based on the original 
Maxwellian distribution were developed afterwards by Pullin fll3 |, Reitz [|Hjj . 



Deshpande and Mandal |D), [IB], [T7[] , Perthame and Coron |TB| , W\, Prendergast 
and Xu [pOR , Xu, Martinelli and Jameson and Raghurama Rao and Deshpande 
22], B3]. For the isentropic Euler system, Kaniel investigated a Kinetic Scheme 



based on an equilibrium distribution function which is different from the classical 
Maxwellian P3L [23]. A general approach to construct equilibrium distributions 



has been presented in [26] and by Perthame 27] who uses an entropy principle. 
For the compressible Navier-Stokes system, Kinetic Schemes were developed by 



Chou and Baganoff 



29, 



A slightly 



and in the group of Deshpande 
different approach was taken by Xu and Prendergast f3lj . 

For scalar conservation laws in one space dimension Backer and Dressier found 
equilibrium distributions following the idea of Kaniel |3^|. In arbitrary space 
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dimensions the scalar case could be treated with a slightly modified transport 
equation [34]]. This approach led to a detailed investigation of the relation 
between the hydrodynamic limit of kinetic equations and nonlinear conservation 
laws by Lions, Perthame and Tadmor [|35|j . 

In this paper, we present a new discrete velocity model based on the methodol- 
ogy of the Kinetic Schemes. In Section |2], the original concept of Kinetic Schemes 
for Euler equations is introduced and then applied with a special equilibrium dis- 
tribution known from the Lattice Boltzmann Method. After that, the obtained 
discrete velocity model is extended to the Navier Stokes case by constructing 
a new discrete Chapman Enskog distribution. In Section [|, consistency of the 
resulting scheme is investigated, leading to an interpretation of both Kinetic 
Schemes and LBM as pseudo-compressibility methods. Section [| concludes with 
numerical results and discussions. 



2 Kinetic schemes for Euler equations 
2.1 Traditional Kinetic Schemes in CFD 

The basis of Kinetic Schemes is the connection between the Boltzmann equation 
of Kinetic Theory of Gases and the macroscopic equations of fluid flow. The fluid 
flow equations can be obtained as (velocity) moments of the Boltzmann equation 

§j- + v-V/ = Q(/). (1) 

Here /(x, v, t) is the velocity distribution function of the gas particles and the 
gradient is taken with respect to the space variable x. The left hand side of the 
equation ([!]) denotes the free flow of the molecules. This free flow is disturbed 
by the molecular collisions, which is represented by the collision term, Q(f), on 
the right hand side of the equation (|]). The mass, momentum and energy of the 
fluid can be obtained as the velocity averages of the particle mass, momentum 
and energy densities. Introducing the notation 

/oo poo poo 
/ / fd Vl dv 2 dv 3 (2) 
-oo J —oo J —oo 

this can be formulated as 

p=(f) pu=(v/> pe=/i|v| 2 A (3) 
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The macroscopic equations can be obtained by integrating the Boltzmann equa- 
tion ([]]), after multiplying it by the vector of the moment functions, as 

^ + v-V/-W)) ) = (4) 

2 i y r / - 

Using (|3|), we get the system 

f t + div (v/> = (Q) 

+ div (v ® v/) = (vQ) (5) 

The mass, momentum and energy are conserved during collisions. Therefore, we 
have 

(Q)=0 (vQ) = ^~|v| 2 Q^ = (6) 

Substituting (|) in @ we obtain 

^ + div (v/> = 

^M + div(v®v/> = (7) 

^M + div/l|v| 2 vA=0 
dt \2 l 1 V 

In the hydrodynamic limit, the gas is dominated by collisions and the particle 
distribution attains the form of a Maxwellian (in the limit of infinite collision 
frequency), given by 

^ ) =(^)T exp (-^) (8) 

This velocity distribution is well known as the one of a gas in (local) thermody- 
namic equilibrium. Hence, the Maxwellian is also called the equilibrium distribu- 
tion. When the distribution is a Maxwellian, the fluxes in (0) can be calculated, 
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yielding 



(v <g> vM) = pu<g>u + pTI and (l;\v\\M \ = p(e + T)u (9) 
where / is the identity matrix. Using the above, we obtain the Euler equations 



as 



at 

d( P e) 
dt 

or equivalently in the form 



^ + div(pu) = 
9(PU) + div(pu ® u + pTI) = (10) 



+ div(p(e + T)u) = 




f = M (11) 



While standard discretizations of the Euler system are based on (|10D, Kinetic 
Schemes use the representation ([□]) which is motivated by Kinetic Theory An 
obvious advantage of (|TTD is the much simpler differential operator which is linear 
and scalar in contrast to the more complicated nonlinear system (|TTj|). 

To discretize (|TTD, traditional CFD techniques like finite difference, finite vol- 
ume, finite element or spectral methods can be applied. Equivalently, one can 
use the Lagrangian approach. In this approach, we replace ( |TTD by the auxiliary 
problem 

|j- + vV/ = 0, f\ t=0 = M (12) 

for which the solution is straightforward, given by 

/(x,v,t) = /(x-vt,v,0) (13) 
Clearly, this solution satisfies 

"as 




a+v-v/iHo. 



However, the constraint / = M. is enforced only initially. The violation of this 
constraint leads to an increasing error as time increases. By stopping the evolution 
after a small time step At and restarting it with a Maxwellian (that has the same 
p, u, e-moments as the solution of the just finished free flow step), the error can 
be kept of order At, giving rise to a first order method for the Euler equations. 
Thus, two clear steps can be identified for the Lagrangian approach: a convection 
step and a relaxation step. In the relaxation step, the velocity distribution relaxes 
completely to the equilibrium distribution. 



2.2 Kinetic Schemes with discrete distributions 

While the Kinetic Schemes mentioned in the above subsection are designed to 
solve the Euler system, it is not necessary to be limited by this restriction. Also, 
the choice of Ai as equilibrium constraint is not mandatory. Obviously, the 
approach is applicable whenever the system of equations allows a representation 
of type (0). In the following, we are going to restrict our considerations to the 
case of isothermal Euler equations, in order to work out the similarities with the 
Lattice Boltzmann method. The isothermal equations (with T = T = (? s where 
c s is the sound speed) are of the form 



or equivalently 



d(pu) 
dt 



1 S \ fdf 



^ + div(pu) = 
+ div(pu <g> u + c 2 s pl) = 



V M ¥ ^-W =0, f=M\ T=To . (15) 



Instead of the classical Maxwellian M. (with fixed temperature T = T ), we 
can also choose any distribution F, as long as the integral expressions in (|15|) 
together with the constraint / = F is equivalent to the Euler system (£[4]). Since 
the integral involves velocity moments up to second order, we are led to the 
following compatibility conditions on the equilibrium distribution 

(F)=P 

(vF) = pu (16) 
(v <g> vF) = pu <g> u + c 2 s pl 
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In particular, we are interested in discrete velocity distributions F which satisfy 
these moment constraints (see also the works of Sanders & Prendergast [[U] 



Nadiga Sz Pullin and Michael Junk (|J7|, PS[). An explicit example in 2D is 
given by the so called D2Q9 distribution used in the Lattice Boltzmann method. 
It is of the form 

8 

F(p,u;v) = ^F J (p,u)(5(v-v J ) (17) 

where 5 is the Dirac-Delta function and 
v = 

v» = V3c s (cos ( (i - 1) f ) , sin ((i- l)f) ) T i = l,... ,4 (18) 
v< = \/6c s (cos ( (i- |) |) , sin ((i - |) f ) ) T i = 5, . . . , 8 

The weights Fi are given by 

F t (p, u) = F*p (l - ^|u| 2 + iu • v, + ^ (u • v,) 2 ) (19) 

with 

F* = - F* = - for % = 1, . . . , 4 = — for i = 5, . . . , 8 (20) 

9 1 9 36 V ; 

In order to obtain a Kinetic Scheme for the isothermal Euler equations, we will 
approximate the equivalent form 

+v ' v/ )> =0 ' f=F (21) 

with the Lagrangian approach described in the last section. Solving the free 
flow equation ^ + v • V x f = 0, starting at time t with equilibrium /(x, v, t) = 
F(p(x, t), u(x, t); v), yields after a time step At 

f(x, v, t + At) = F(p(x - vAt, t), u(x - vAt, t); v) 

8 

= F M* - vAt, t), u(x - vAt, t))5(v - v,). 

i=0 
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Using the relation ip(y)S(y — Vj) = i/j(vi)5(v — Vj), which holds for any continuous 
function ■0, we obtain further 

8 

/(x, v, t + At) = F i(p(* ~ v i At ' *)» u ( x " V * A ^ *))<K V - v *)- 

i=0 

Denoting the weights of the discrete distribution /(x, v, t + At) by /j(x, t + At), 
the evolution can also be described without mentioning the Dirac deltas at the 
(fixed) discrete velocities 

fi(x,t + At) = F(p(x-v i At,t),u(x-v J At,t)), i = 0,...,8. (22) 

Since integer multiples of v^At make up a regular square grid which is invariant 
under VjAt-translations, the scheme ( |2"2"D only accesses nodal data if x is also a 
node of the grid. We remark that the grid length is given by 

Ax = V3c s At. (23) 



The connection between fl22|) and the classical Lattice Boltzmann method be- 
comes most obvious under the change of variables x i— ► x + VjAt, which leads 
to 

/ i (x + v l At,t + At) = F i (p(x,t),u(x,t)), i = 0,...,8. (24) 
Indeed, ([23]) coincides with the Lattice Boltzmann evolution [^, ^ 



fi (x + v, At, t + At) - fi (x, t) = — ^ (F (x, t) - /< (x, t)) (25) 

t« 



if we set t^ = At (see references |39| and for the LBM based on BGK- model). 



At first glance, this seems to be a contradiction, because the Kinetic Scheme has 
been set up for the Euler system while it is known that the Lattice Boltzmann 
method approximates the Navier Stokes system. In fact, setting = At amounts 
to a high dose of viscosity (typically, LBM applications are run with £r in-between 
| At and At). The apparent contradiction is resolved with the remark that the 
Lagrangian approach to Kinetic Schemes yields only a first order method. The 
numerical viscosity of that scheme is quite high, particularly in applications with 
low Mach number flows. This numerical viscosity of the Kinetic Scheme is exactly 
the viscosity corresponding to £# = At in LBM and thus has a physically correct 
structure. In |38|| , a Kinetic Scheme for the Euler system could therefore be used 



as solver for the Navier Stokes equations. Huang et al. [[IT], in their Lattice 



Boltzmann Method for compressible flows, used a similar approach, even though 
they did not mention Kinetic Schemes. 
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3 Extension to Navier— Stokes equations 



Kinetic Schemes can also be extended to the case of Navier-Stokes equations, by 
using a Chapman-Enskog distribution function Tce instead of the Maxwellian 
constraint M. in (|i"4|). This approach has been pursued in |29], The 
Chapman-Enskog distribution function Tce is obtained as a small perturbation 
of the Maxwellian. See references [i2, 53, 15] and [23] for details of the derivation. 
For a mono-atomic gas in three dimensions, the distribution function is of the 
form 



T 



CE 



where 



-K 



M 



dT 

dxi 



p 2T° tCj p T Ci 



5 2T J 



P 



dui On, 



2 . duk 



dxj dx 



3 6ij dx k 



(26) 



(27) 



and c = v — u is the peculiar velocity. Here, we will again consider the simpler 
case of isothermal equations in 2D. Following [151 
form 



the equations are of the 



^ + div(pu) 



0. 



d(pu) 
dt 



(2* 



+ div(pu ® u + c 2 s pl) = div r\ 



where the viscous stress tensor is given by 
r] = vp{2S + div ul), 



duj dm 



+ 



dxj dxj 



(Note that div 7] is the vector obtained by applying divergence to the rows of rj.) 
Equivalently, we can write system (EHI) as 



df 
dt 



+ v-VJ 



0. 



/ 



CE 



(29) 



where Fqe satisfies the moment constraints 

(F C e) = P 
(vF CE ) = pu 
(v ® vF C e) = pu <g> u + c 2 s pl 



(30) 
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Instead of using the continuous distribution function Tqe (given in (^) with 
qi = and T = T ), we construct a discrete F C e satisfying (|30|). One possibility 
is to discretize the classical Chapman-Enskog distribution function in the velocity 
variable. A Kinetic Scheme based on this approach will be presented elsewhere 



46fl . Here, we follow a different idea based on a general solution technique for 



moment problems of the form (|30|) which uses orthogonal polynomials | 26| . For 
the special D2Q9 model, however, it is not necessary to work out the general 
ideas. In fact, conditions (|30"D can be reduced to those of the Euler system if we 
replace puCg>u by pu®u — 77. This observation can be used, if we write the weights 
( |T9| ) of the equilibrium distribution function ([H]) in terms of pu <g) u. Introducing 
the matrix product 

2 



we find 



2 

2 



and 



u ® u : v ® v = UiUjViVj = (u ■ v 

2 

u ® u : / = UiUjSij = |u| 2 



Fi(p, u) = F*p ( 1 + — u • v< + <g> u : ( — v< ® v< - I 



so that Ql9|) can be written as 

1 1 / 1 

-u.v l + -u^u:^- 

Replacing u®uby u®u - z/(2S' + div uJ), we finally obtain 

^C£,i(P,u) 

= i^*p ( 1 + -^u ■ v, + ^(u ®u-2vS - z/divu/) : ( — v< ® v, - / ) ) . (31) 



,s 



or after going back to scalar products in and u, 

F C E,i{p, u) = i^P^l + • Vi - ^|u| 2 + ^i(u • Vi) 2 

v „ 2/ -. / .l , |2 



S:V(0Vj _ divu | Vj |^_ 2 . (32 ) 
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It is easy to check that the so defined Chapman-Enskog distribution 

8 

F CE (p, u; v) = ^ F CE,i{pi u )^( v - v *) ( 33 ) 

i=0 

satisfies fl3~0|) . We also remark that F C e is a perturbation of the original D2Q9 
equilibrium distribution, similar to the classical case in Kinetic Theory where 
fl26p is a perturbation of the Maxwellian (||). 

To develop a Kinetic Scheme for Navier-Stokes equations, we follow the same 
procedure as in the previous section, except that the distribution used as con- 
straint after every time-step will now be the Chapman-Enskog distribution, Fqe- 
We end up with the scheme 

/<(x,t + At) =F CE ,i{p{x- v;AM),u(x- ViAt,t)), i = 0,..., 8 (34) 

where the moments are updated according to 

8 

p(x,t + At) = ^2fi(x,t + At) 

1=0 „ (35) 
u(x,t + At)= / + gv^t + At). 



4 The incompressible limit 

To investigate the behavior of the Kinetic Scheme at low Mach numbers, we first 
scale the compressible Navier Stokes system appropriately Low Mach number 
flows appear if u is very small compared to c s . Taking a typical speed U and 
length scale L of the flow, the time scale G is chosen in accordance to these scales 
as 




The density p is assumed to be of order one so that no scaling is needed. To avoid 
superscripts, we will not change the symbols for scaled functions and variables. 
If we refer to unsealed quantities (which appear less often in this section), we add 
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a hat to the symbols. After some algebra, we obtain the scaled version of 

- + — div(pu) = 0, 
d(pu) QU . . c 2 9 

By assumption, QU/L = 1 and c 2 Q/(LU) = c 2 /U 2 . Introducing the Mach num- 
ber Ma = U/c a and the Reynolds number Re = UL/v of the flow, we end up 
with 

dp 

+ div(pu) = 0, 

d(pu) 1 1 ( 36 ) 

+ div(pu ® u) + — Vp = — div(2p5 + p div uJ) 

If we approximate (|36|) by the Kinetic Scheme introduced in then the time 
step is related to the grid length by At = V3c s Ax (see or in scaled quantities 

At = v^MaAx. (37) 

This relation already indicates the typical problem that any explicit solver for 
the compressible equations faces in the incompressible limit: to get a reasonable 
space resolution, the time resolution must be extremely fine (if Ma <C 1) to satisfy 
the CFL-condition (|37|). 

To find out which equations are approximated by the Kinetic Scheme, we 
perform a consistency analysis in the coupled limit At, Ma — > 0. More precisely, 
we assume 

~' A = const for At -> 0, Ma -> 0. (38) 



Ma 2 



To begin with, let us rewrite the Chapman Enskog distribution (|3lD in scaled 
quantities. 

F CE ,i(p, u) = F*p ^1 + Ma 2 u • v< 

Ma 2 / 1 \ 

+ _U®„-_(2S + divuO) : (Ma 2 Vi®Vi-/) |. (39) 
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Setting F C e{p, u ; v ) = Yli=o FcE,i{p, u)5(v — Vj) and using p(x) and u(x) as 
initial values, the Kinetic Scheme yields at the end of the first time step 

^(x) = (F C£ (p(x - vAt), u(x - vAt); v)) (40) 

and 

(p 1 u 1 )(x) = (vF os (p(x-vAt),u(x-vA*);v)). (41) 

To obtain a Taylor expansion around At = we need At-derivatives of ( fffl| ) 
and ( pfTD up to a certain order. Obviously, each At-derivative leads to a space 
derivative with — v as factor (i.e. j^i = ~ Vi '£r)- ^° ^ rs ^ or der consistency 
in At, we nevertheless need higher At-derivatives. This is due to the fact that 
terms of the form At 2 /Ma 2 , At 2 /Ma 3 and At 3 /Ma 4 are not negligible in the 
coupled limit (|3^). Consequently, we also need higher order v-moments of the 
Chapman Enskog distribution. Taking the scaling into account, we get from (3C) 

(F C e) = P 

(vF CE ) = pu ( 42 ) 
(v <g> wF CE ) = pu ® u + ^pl - ^(2p5 + pdivu/). 

The third order moment can be calculated using the explicit form of Fqe given 
in (§§. We find 

(viVjV k F C E) = ^p^p(6ijUk + ^ ikU i + &kjUi). (43) 

Finally, from the fourth and fifth order moments we only need to know the terms 
of leading order 



(viVjVkViFcE) = + tikSji + 8u8jk) + O y^^J 

(ViVjVkViVmFcE) = O 

The Taylor expansion of (HU1) is then given by 



Ma 4 



(44) 



P 1 = (Fas) - £ (v iFcE ) At + \^*- (viVjFcE) At 2 



1 d 3 

i a a a (viVjV k FcE) At 3 + 

O OXiOXjOXk 
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Since 



At 2 

(v^Fce) At 2 = ^pSij + 0{At 2 ) 
, , n ( At 3 \ , n , . , 



v l v j v k F CE ) At 6 = O [y^j = 0(At 2 ) (45) 



we conclude 



(viVjvmFcE) At 4 = O (j^) = 0{At 2 ) 



p 1 = p + ( -XAp - div(pu) ) At + 0{At 2 ). (46) 



Similarly, we get for the momentum defined in (fy|) 
(pV), = ( Vi F C e) - 4- (v tVl F CE ) At + I— — - (wjVtFcB) At 2 



dxi 2 dxidxj 



d 3 

(viVjVkViFcE) At 3 + . . . 



6 dxidxjdxk 

While the second order moments yield exactly the fluxes of momentum, the third 
order moments give rise to some additional terms. Using (ff3|), 



According to fl44]), the fourth order moment leads to 

1 d 3 , „ v A 3 1 d A A , , m ( At 3 



(v^v^Fce) At d = --— ApAt + O 



6dxidxjdxk 2 dxi \Ma 2 

and fifth order moments are negligible since 



(ViVjVmVmFoE) At 4 = O [ ) = 0(At z 



Thus 



P l u] = pUl + \^XA{pui) + A^- div(pu) - ^^Ap 



(47) 
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Since we assume that all appearing quantities are scaled, the equation for p 1 u 1 can 
only be balanced if dp/dxi = C(Ma 2 ). Hence, we assume that p = p(l + Ma 2 p) for 
some constant p > and a function p which is assumed to be of order one together 
with its derivatives in the limit under consideration. (This is the standard scaling 
of the density in isothermal, low Mach number flows.) Using the additional 
information on p and observing that Ma 2 = O(At), we can simplify (f46|) in 
lowest order to 

divu = C(At). (48) 

This equation has to be understood in the sense that the order one assumption 
on u and p is only consistent if the divergence of u is O(At). Before we explain 
how the Kinetic Scheme guarantees the approximate divergence-free condition, 
we note that relation (^) and the structure of p reduces (||7|) to the Navier Stokes 



equation with a first order error term 



^ + (u.V)u + Vp=(^- e + h)Au+ O(At). (49) 



To explain the mechanism that leads to ([48"D, we use (^) again, keeping the first 
order terms. After division by At and Ma 2 this leads to 

d P v / x 1 v 1 . . At 2 s 

+ div(pu) + — — div u = -XAp + 0(- 



dt vr y Ma 2 2 ^ v AtMa 2 '' 

To resolve the additional terms of order one on the right hand side, we have 
to expand (H5[) one order higher. Using the explicit knowledge of the relevant 
moments, relation (48) and our assumption on p, we find 



dr> 11 

+ (u ■ V)p + ^2 div u = -XAp + div ((u ■ V)u) + 0{At). 



Note that equations of this type are used in pseudo-compressibility methods E7 



[48[ to ensure the divergence free condition. In fact, it uses elements of Chorin's 
artificial compressibility method to replace div u by the equation 

dp 

e— + div u = 

at 

and of the pressure stabilization method 

div u — e Ap = 
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which was originally used by Hughes, et. al., [pD| . However, the convection term 
and the nonlinear term which follow automatically from the kinetic approach are 
usually not considered. We thus conclude that: 

In the coupled limit At, Ma — > with At/ Ma 2 = A with the assump- 
tion that p = p(l + Ma 2 p) and that u, p and their derivatives are order 
one functions, the Kinetic Scheme is consistent to the incompressible 
Navier Stokes equation with effective Reynolds number 

J_ J_ b. 

~Re! ~ Ye + 2' 

The scheme can be viewed as a new pseudo-compressibility method. 

Note that in the case Re = oo, the Chapman Enskog distribution reduces to 
a Maxwellian and the Kinetic Scheme is equivalent to the Lattice Boltzmann 
method with relaxation parameter tn = At. Therefore, LBM can also be viewed 
as a pseudo-compressibility method in that case. Since an additional viscosity 
term appears in the coupled limit At, Ma — ► 0, the Kinetic Scheme with v — still 
approximates the solution of an incompressible Navier Stokes equation. As al- 
ready mentioned earlier, this idea has been used in [|38| to construct Navier Stokes 
solutions with a Kinetic Scheme which is just based on a discrete Maxwellian. 

5 Numerical Results and Discussions 

We first note that the term involving div u in the Chapman Enskog distribution 
( |32| ) is actually not important in low Mach number situations and thus can be 
neglected. Note that such modifications are very simple in the framework of 
Kinetic Schemes: by adding or deleting terms in the distribution function, the 
macroscopic equations can easily be modified. In the case of LBM, on the other 
hand, the Chapman Enskog distribution is implicitly given through properties of 
the collision operator which makes it much harder to develop such schemes for 
modifications of the incompressible Navier Stokes equation. 

We adjust the viscosity parameter v in the Chapman Enskog such that the 
effective viscosity turns into the required one. This prevents the numerical viscos- 
ity from spoiling the results of the simulations. Altogether, we base our Kinetic 
Scheme on the following distribution function (which is now written again in 
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unsealed variables) 



U • V, 1 , ,n 1 - ,n / Z/ At\ „ V 



In a first test case, we apply the scheme to a Poiseuille flow in an infinitely 
long channel (in a^-direction) of width one with a constant acceleration g. The 
incompressible Navier Stokes solution for this case is explicitly known to be 

ui(x 2 ) = — (1 - x 2 )x 2 , u 2 = (50) 
2v 

with a constant pressure. In our simulation we choose p = 1 initially. The in- 
finitely long channel is modeled by periodic boundary conditions in xi-direction. 
The fixed wall conditions for u are enforced simply by setting u = at the bound- 
ary nodes. In contrast to LBM, where the no-slip condition has to be enforced 
by properly setting the incoming occupation numbers, no such complications are 
found here, because the unknowns in the Kinetic Scheme are directly the flow 
variables p and u. The boundary conditions for density can be obtained from 
the Navier Stokes equation (|28|). Multiplying the equation with the outward unit 
normal vector n and observing that u = at the boundary, we get the condition 

(? S -Jr = n • divr]. (51) 
on 

For the exact solution (|5C|) one easily checks that 

divr? = -vpg 

so that n • div r) — at the upper and lower walls giving rise to homogeneous 



boundary conditions for p. (According to |47| , homogeneous Neumann conditions 
for p are also reasonable in more general, moderate flow situations.) The force 
term is incorporated into our scheme by a splitting approach: in a first step the 
Kinetic Scheme approximates the Navier Stokes evolution and in a second step, 
the acceleration is taken care of by an explicit Euler step for the velocity variable. 
To calculate the stress tensor S, we use central differences. 

^From the solution (0), we can see that the maximum velocity 

u = ± 

8u 
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is obtained at the center of the channel. By setting g = 0.01 and v = 0.01, we 
get U = 0.125 (note that U is the Mach number since c s = 1). With 11 points 
across the channel and initial velocity u = 0, we find a numerical approximation 
which reproduces the predicted parabolic shape (see Fig. [l|). The other velocity 
component stays zero and the density remains constant. 
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Figure 1: Poiseuille velocity profile 

Due to the symmetries in this simple test case, the incompressibility condition 
is satisfied exactly. Consequently, compressibility errors are not present and the 
accuracy of the scheme just depends on how closely the steady state is approxi- 
mated. For several values of t, the L°° error behaves as depicted in Fig. 0. In all 
calculations, the number of grid points is 11 and U = 0.125. 
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Figure 2: L°°-error versus time 

The next test case is a slight modification of the previous one, where the top wall 
of the channel now moves with a fixed velocity w in ^-direction (Couette flow). 
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Here, the exact solution differs from ( |5*0"D by an additional linear term 



9 



ui(x 2 ) = — (1 - x 2 )x 2 + wx 2 , 



u 2 = 0. 



(52) 



Again, the u-boundary condition at the moving wall is easily enforced by setting 
ui = w and u 2 = 0. Using the same settings as above with w = 0.12 the results 
are again in agreement with the exact solution (see Fig. |3|). 
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Figure 3: Velocity profile for Couette flow 

Our final test case is the driven cavity flow problem. The incompressible fluid 
is now bounded by a square enclosure with side length one. The top lid, which 
moves with velocity U, generates the fluid motion in the cavity which shows 
typical vortex phenomena. For our calculations we use a 129 x 129 uniform grid 
and Reynolds numbers 100, 400 and 1000. The lid velocity U is set to one and 
c s = 10 in all cases. The calculations are initialized with p — 1 and u = inside 
the cavity. As termination criterion we choose a residue fall of 3.75 decades in 
the equation for p. The typical number of cycles to get steady state solutions 
is 100,000. We remark that no special attention has been paid to acceleration 
of convergence. Our aim is only to show that the new discrete velocity method 
works for complex test problems. (Note that the pressure develops singularities 
in the top corners due to the jump in the boundary conditions for velocity.) 
In order to demonstrate that the Kinetic Scheme can be used like a Lattice 
Boltzmann method we implement the boundary conditions using the fast bounce 
back algorithm f51 |. To explain this approach we remark that on the kinetic level, 



boundary conditions are required for the transport part of the equation 

df 
dt 



+ v.V/ = 0. 



(53) 
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Since fl53|) is a linear hyperbolic equation, information has to be provided for 
those characteristics which enter the domain at a boundary. In our model, the 
characteristics are straight lines along the discrete velocity directions vi...v§. 
The bounce back condition sets the value for the information of an incoming 
direction equal to the information that leaves the domain in the opposite direc- 
tion, which is easily available due to the symmetry of the discrete velocity set. It 
can be shown |51] that these conditions simulate no slip conditions at the Navier 
Stokes level. At the upper lid, a modification is required which takes care of the 
momentum flux generated by the movement fl52| |. To illustrate our results, the 
horizontal velocity component Mi is shown along a vertical section through the 
center of the cavity (Fig. ^). Similarly, we plot the vertical component u<i along 
the central horizontal section (Fig. The results are compared with those ob- 
tained by Ghia et. al. [53 and they are in good agreement. In Figures ^ and ||, 
the symbols refer to the tabulated simulation results in and the lines refer to 
the results obtained by the new Kinetic Scheme. 
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Figure 4: wi-velocity along a vertical line through the center of the cavity 
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Figure 5: ^-velocity along a horizontal line through the center of the cavity 

Plots of the stream functions are given in Figures ^ to [|. We remark that the 
stream function ip is, strictly speaking, not well defined because the approximate 
velocity field is not exactly divergence free. In this problem is discussed for 
the Lattice Boltzmann method and we use the proposed numerical procedure for 
the calculation of ip (integration of U2 along horizontal sections from left to right). 
The levels of the isolines are those from ||53|1 . We limit ourselves in this study to 
the use of uniform grids, as our purpose is to show that the new discrete velocity 



model works. With clustered grids, the solution can be different 54 . 
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Figure 8: Re = 1000 



The numerical cost of the Kinetic Scheme is directly comparable to that of the 
Lattice Boltzmann method based on the D2Q9 model. Both algorithms have the 
same structure consisting of a propagation and a collision step. The only differ- 
ence is that in the Kinetic Scheme the stress tensor has to be calculated (by taking 
central differences of the velocity field) and that the equilibrium distribution is 
extended by the viscosity term. On the other hand, the Kinetic Scheme needs 
less memory because there is no need to store the occupation numbers. Apart 
from two copies of p and u (new and old time step) an efficient implementation 
requires three more variables to store the stress tensor. Altogether a 2D compu- 
tation needs nine floating point variables per node, independent of the underlying 
discrete velocity model. Compared to that, D2Q9 Lattice Boltzmann methods 
need 21 variables per node (for p, u and two copies of the occupation numbers 
/o, . . . , fs) and the number increases if models with more velocities are used. Also, 
when passing over to 3D calculations, the memory usage of the Kinetic scheme 
increases by five variables per site whereas D3Q15 Lattice Boltzmann methods 
need 13 more variables in each node. Of course, the discrepancy becomes even 
larger if multi-phase flows are simulated. Then, for simple algorithms, the mem- 
ory requirement has essentially to be multiplied by the number of participating 
species. Taking these considerations into account, Kinetic Schemes seem to be a 
powerful alternative to Lattice Boltzmann methods. On the one hand, they are 
formulated in the same kinetic framework allowing the use of LBM specific solu- 
tion techniques (kinetic boundary conditions, treatment of phase boundaries in 
multi-phase flows, etc.). On the other hand, Kinetic Schemes only use the actual 
flow variables and thus can profit directly from established Finite Difference or 
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Finite Volume methods. In addition, the memory consumption is greatly reduced 
compared to Lattice Boltzmann algorithms. 

Since the idea of LBM is to use a kinetic model which is as simple as pos- 
sible under the constraint that the macroscopic limit equations are correct, the 
method is not capable of quantitatively predicting the behavior of a rarefied gas 
and should therefore only be applied close to equilibrium situations. To show 
consistency of LBM to the Navier Stokes equation, exactly this equilibrium as- 
sumption is used in the Chapman Enskog expansion which amounts to assuming 
that the occupation numbers are given by a Chapman Enskog distribution. A 
natural idea is therefore to build the Chapman Enskog distribution directly into 
the algorithm which is exactly the construction principle of the present Kinetic 
Scheme. Thus, Kinetic Schemes can be viewed as a consequent advancement of 
the Lattice Boltzmann Method. 

We conclude our discussion with a remark concerning the extension to the full 
Navier Stokes system including the energy equation. A fundamental problem of 
the basic Lattice Boltzmann method based on a simple BGK collision operator is 
that the Prandtl number is not a free parameter. In a Kinetic Scheme, the heat 
conduction and viscosity parameters enter directly into the Chapman Enskog 
distribution (similar to the continuous case fl27D ) and thus can naturally be varied 
independently. 

6 Conclusions 

The similarities and differences between the Lattice Boltzmann Method, which 
has recently become popular, and the Kinetic Schemes, which are routinely used 
in Computational Fluid Dynamics, are studied. A new discrete velocity model for 
the numerical simulation of incompressible Navier-Stokes equations is presented 
by combining both the approaches. This approach of Kinetic Schemes with dis- 
crete distributions is shown to be more convenient and useful compared to the 
Lattice Boltzmann Method. Since both methods coincide for a particular choice 
of parameters, the analysis of the Kinetic Scheme also applies directly to LBM 
in that case. In particular, the conclusion that the Kinetic Scheme is a special 
pseudo-compressibility method illuminates the Lattice Boltzmann approach. 
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